
library(Hmisc)
library(ggplot2)
library(foreign)

setwd("") #working directory with R code
stem <- "" # path for data files


multiplot <- function(...,title=title,plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)
  
  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)
  
  numPlots = length(plots)
  
  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols) + cols * 1),
                     ncol = cols, nrow = (ceiling(numPlots/cols) + 1))
  }
  
  if (numPlots==1) {
    print(plots[[1]])
    
  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout), heights = unit(c(1,4,4,4), "null"))))
    grid.text(title, vp = viewport(layout.pos.row = 1, layout.pos.col = 1:3))
    # Make each plot, in the correct location
    
    index <- c(2,3,4,6,7,8,10,11,12)
    
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == index[i], arr.ind = TRUE))
      
      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}

plotting <- function(data,title,s_vtc,s_cddc){

  vals <- 1:365
  
  # Initialize lists
  plots <- vector("list", 9) 
  plots2 <- vector("list", 9) 
  diff <- vector("list", 9) 
  
  for (i in 1:9) {
    mydata <- subset(data,data$combin==i)
    mydata_vtc <- subset(mydata,mydata$site==0)
    mydata_cddc <- subset(mydata,mydata$site==1)
    
    # span = 0.4 best for CDDC, 0.6 best for VTC
    loess_vtc <- loess(ls2 ~ x, span=s_vtc[i], data=mydata_vtc)
    loess_cddc <- loess(ls2 ~ x, span=s_cddc[i], data=mydata_cddc)
    
    predict_vtc <- predict(loess_vtc, data.frame(x=vals))
    predict_cddc <- predict(loess_cddc, data.frame(x=vals))
    
    temp3 <- predict_vtc - predict_cddc
    diff[[i]] <- data.frame(x = vals , y = temp3 , combin = i)
      
    temp <- cbind(y =  predict_vtc , site = 0 , x=vals)
    temp2 <- cbind(y =  predict_cddc, site = 1, x=vals)
    predict <- as.data.frame(rbind(temp,temp2))
    
    plots[[i]] <- ggplot(data = mydata, aes(y = ls2, x = x,colour=factor(site))) + 
                  geom_point() +
                  geom_line(data = predict, aes(y = y, x= x,colour=factor(site))) + 
                  theme(legend.position="none") + 
                  ggtitle(paste(titles[[i]])) +
                  labs(x="Time since release/discharge (days)",y="Log(-Log(Survival))") + 
                  theme(plot.title = element_text(family = "Times", color="#666666", face="bold", size=10, hjust=0)) +
                  theme(axis.title = element_text(family = "Times", color="#666666", size=10)) 
    plots2[[i]] <- ggplot(data = diff[[i]], aes(y = y, x= x)) + 
                    geom_point() 
  }
  
  # pdf("~/plots.pdf")
  
  multiplot(plots[[1]], 
            plots[[2]], 
            plots[[3]], 
            plots[[4]], 
            plots[[5]], 
            plots[[6]], 
            plots[[7]], 
            plots[[8]], 
            plots[[9]],title=title,cols=3)
  
  # dev.off()
  
  multiplot(plots2[[1]], 
            plots2[[2]], 
            plots2[[3]], 
            plots2[[4]], 
            plots2[[5]], 
            plots2[[6]], 
            plots2[[7]], 
            plots2[[8]], 
            plots2[[9]],title=title,cols=3)
  
  diffs <- rbind(diff[[1]],
                diff[[2]],
                diff[[3]],
                diff[[4]],
                diff[[5]],
                diff[[6]],
                diff[[7]],
                diff[[8]],
                diff[[9]])
  
  return(diffs)

}

s_vtc <- c(rep(.6,36),rep(.7,9),rep(.6,9),rep(.7,9),rep(.5,9),rep(.7,9),rep(.7,9),rep(.5,9),rep(.7,9),rep(.7,27))
s_cddc <- c(rep(.4,36),rep(.7,36),rep(.6,36),rep(.9,9),rep(.8,9),rep(.8,9))

titles <- vector("list", 9) 

titles[[1]] <- "CDDC(SR) + VTC (SR)"
titles[[2]] <- "CDDC(SV14) + VTC(SR)"
titles[[3]] <- "CDDC(SV0) + VTC(SR)"
titles[[4]] <- "CDDC(SR) + VTC(SV90)"
titles[[5]] <- "CDDC(SV14) + VTC(SV90)"
titles[[6]] <- "CDDC(SV0) + VTC(SV90)"
titles[[7]] <- "CDDC(SR) + VTC(SV30)"
titles[[8]] <- "CDDC(SV14) + VTC(SV30)"
titles[[9]] <- "CDDC(SV0) + VTC(SV30)"

var <- c(rep("opi",4),rep("amp",4),rep("any",4),rep("mtd",3))
type <- c(rep(c("A","B","C","H"),3),c("A","B","C"))
temp <- vector("list", 15) 

for (i in 1:15) {
  
  fname <- paste(stem,"mydataY",type[i],"_",var[i],".xpt",sep="")
  data <- sasxport.get(fname)
  
  title <- paste("Y",type[i],"_",var[i],sep="")
    
  colnames(data)[2] <- "x"
  colnames(data)[1] <- "site"
  
  temp[[i]] <- plotting(data,
                        title,
                        s_vtc[((i-1)*9+1):(i*9)],
                        s_cddc[((i-1)*9+1):(i*9)])
  
  temp[[i]]$var <- var[i]
  temp[[i]]$type <- type[i]
  
}

new_data <- rbind(temp[[1]],
                  temp[[2]],
                  temp[[3]],
                  temp[[4]],
                  temp[[5]],
                  temp[[6]],
                  temp[[7]],
                  temp[[8]],
                  temp[[9]],
                  temp[[10]],
                  temp[[11]],
                  temp[[12]],
                  temp[[13]],
                  temp[[14]],
                  temp[[15]])

fname2 <- paste(stem,"diffs.dbf",sep="")

write.dbf(new_data,fname2)
